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Abstract. This paper describes the variational formulation of tem- 
plate matching problems of computational anatomy (CA); introduces 
the EPDiff evolution equation in the context of an analogy between CA 
and fluid dynamics; discusses the singular solutions for the EPDiff equa- 
tion and explains why these singular solutions exist (singular momentum 
map). Then it draws the consequences of EPDiff for outline matching 
problem in CA and gives numerical examples. 

"I shall speak of things ... so singular in their oddity as in some manner 
to instruct, or at least entertain, without wearying." - Lorenzo da Ponte 



1 Introduction 

Computational Anatomy (CA) must measure and analyze a range of variations 
in shape, or appearance, of highly dcformablc structures. The problem statement 
for CA was formulated long ago 

In a very large part of morphology, our essential task lies in the compar- 
ison of related forms rather than in the precise definition of each. . . . 
This process of comparison, of recognizing in one form a definite permu- 
tation or deformation of another, . . . lies within the immediate province 
of mathematics and finds its solution in . . . the Theory of Transforma- 
tions. - D'Arcy Thompson, On Growth and Form (1917) 

The pioneering work of Bookstein, Grenander and Bajscy |2I3I4| first took up 
this challenge by introducing a method called template matching. The past sev- 
eral years have seen an explosion in the use of template matching methods 
in computer vision and medical imaging that is fulfilling D'Arcy Thompson's 
expectation |RK{nriK)lll)lllliai:^ll4llRll<{ll7llrill^ . These methods enable the 
systematic measurement and comparison of anatomical shapes and structures in 
medical imagery The mathematical theory of Grcnander's dcformablc template 
models, when applied to these problems, involves smooth invertible maps (dif- 
feomorphisms) , as presented in this context in |9I1UI18I19I2UI21| . In particular, 



the template matching approach involves Riemannian metrics on the diffeomor- 
phism group and employs their projections onto specific landmark shapes, or 
image spaces. 

The problem for CA then becomes to minimize the distance between two 
images as specified in a certain representation space, V. Metrics are written 
so that the optimal path in V satisfies an evolution equation, which was first 
discovered in abstract form [22| and later called EPDiff when it arose in the 
Euler-Poincare theory of optimal motion on smooth invcrtible mappings called 
diffcomorphisms, |23| . 

The EPDiff equation coincides with the Eulcr equation for ideal fluids in the 
case that the Riemannian metric for the distance between two images is the L 2 
norm. Another type of norm on V (called the H 1 norm) arises in the theory of the 
fascinating nonlinear coherent solutions of shallow water waves called solitons. 
Solitons interact with each other elastically, so they re-emerge unscathed from 
fully nonlinear collisions. EPDiff with the H 1 norm on V describes the peaked 
soliton solutions of the Camassa-Holm shallow water wave equation. As we shall 
see, the Camassa-Holm peakons arise from a general property of Hamiltonian 
systems called their momentum map. A discussion of EPDiff and peakons in 
the particular case of template matching appears in |24| . 

In this paper, we shall draw parallels between the two endeavors of fluid 
dynamics and template matching for computational anatomy, by showing how 
the Euler-Poincare theory of ideal fluids can be used to develop new perspectives 
in CA. In particular, we discover that CA may be informed by the concept of 
weak solutions, solitons and momentum maps for geodesic flows |24I25I26| . 

1.1 Problem &: Approach for Computational Anatomy 

Computational Anatomy (CA) compares shapes (graphical structures) by mak- 
ing a geodesic deformation from one shape to the another. Among these 
graphical structures, landmarks and image outlines in CA are found to be sin- 
gular solutions of the geodesic EPDiff equation. A momentum map for sin- 
gular solutions of EPDiff yields their canonical Hamiltonian formulation, which 
provides a complete parameterization of the landmarks and image outlines 
by their canonical positions and momenta. The momentum map provides 
an isomorphism between landmarks (and outlines) for images and singular 
(weak) solutions of EPDiff. (These solutions are solitons in ID.) This isomor- 
phism provides for CA: (1) a complete and non-redundant data representa- 
tion; (2) a dynamical paradigm in which image outlines interact by exchange 
of momentum; (3) methods for numerical simulation & data assimilation. Euler- 
Poincare theory also provides a framework for unifying and extending the various 
approaches in CA. 

Thus, the concept of momentum becomes important for CA, because mo- 
mentum: 

— Completes the representation of images (momentum of cartoons); 

— Informs template matching of the possibility of soliton-likc collisions and mo- 
mentum exchange in image outline interactions; 



— Encodes the subsequent deformation into the initial locus and momentum 
of an image outline; 

— Provides numerical simulation methods using the momentum map for 
right action as a data structure; and 

— Accomplishes matching and data assimilation via the adjoint linear prob- 
lem for template matching, using the initial momentum as a control vari- 
able. 

All of these momentum properties flow from the EPDiff equation. 

Outline of the paper. Section [21 describes the template matching variational 
problems of computational anatomy, explains the analogy with fluid dynamics 
and introduces the fundamental EPDiff evolution equation. The singular solu- 
tions for the EPDiff equation 12.1fl with diffcomorphism group G are discussed 
in section|3| They are, in particular, related to the outline matching problem in 
computer vision, examples of which are given in section 

2 Mathematical formulation of template matching for CA 
2.1 Cost 

Most problems in CA can be formulated as: Find the deformation path 
(flow) with minimal cost, under the constraint that it carries the tem- 
plate to the target. Such problems have a remarkable analogy with fluid dy- 
namics. The cost assigned in template matching for comparing images Xq & X\ 
is specified as a functional 



defined on curves ipt in a Lie group with tangents = Ut o cp t and 2 t = tpt • 2o- 
In what follows, the function u t ^ £(u t ) = ||u t || 2 will be taken as a squared 
functional norm on the space of velocity vectors along the flow. The Lie group 
property specifies the representation space for template matching as a manifold 
of smooth mappings, which may be differentiated, composed and inverted. The 
vector space of right invariant instantaneous velocities, u t = {dft/dt) o ip^ 1 
forms the tangent space at the identity of the considered Lie group, and may be 
identified as the group's Lie algebra, denoted g. 

2.2 Mathematical analogy between template matching and fluid 
dynamics 

(I) The frameworks of CA and fluid dynamics both involve a right-invariant 
stationary principle with action, or cost function. The main differences arc that 
template matching is formulated as an optimal control problem whose cost func- 
tion is designed for the application, while fluid dynamics is formulated as an 
initial value problem whose cost function is the fluid's kinetic energy. 




(II) The geodesic evolution for both template matching and fluid dynamics is 
governed by the EPDiff equation |27I21| . 

(J^ + u • V)m+ (Vu) T -m + m(divu) =0. (2.1) 

Here u = G * m, where G* denotes convolution with the Green's kernel G for 
the operator Q op , where 

m = — =: Q op u 
ou 

The operator Q op is symmetric and positive definite for the cost defined by 

Cost(i h-> (p t ) = ( t{u t )dt=\ ( \\vi t \\ 2 dt = \ [ (u t ,Q op u t )dt 
Jo z Jo ^ Jo 

with L 2 pairing ( • , • ) whenever ||u t |j 2 is a norm. 

(III) The flows in CA and fluid dynamics both evolve under a left group action 
on a linear representation space, T t = (ft ■ 1o- They differ in the roles of their 
advected quantities, at — clqo ip^ 1 . The main difference is that image properties 
are passive and affect the template matching as a constraint in the cost function, 
while advected quantities may affect fluid flows directly, for example through the 
pressure, so as to produce waves. 

2.3 How EPDiff emerges in CA 

Choose the cost function for continuously morphing Xq into I\ as 

Cost(i h+ (p t ) = [ £(vL t )dt= f \\u t \\ 2 dt, 
Jo Jo 

where it t is the velocity of the fluid deformation at time t and 

llutll 2 = (u t , QopUt } , 

and Q op is our positive symmetric linear operator. Then, the momentum gov- 
erning the process, m t = Q D pU t with Green's function G : u t = G * m t satisfies 
the EPDiff equation, (j2.1(l . This equation arises in both template matching and 
fluid dynamics, and it informs both fields of endeavor. 

2.4 Deriving EPDiff from Euler-Poincare Reduction of Hamilton's 
principle 

Euler-Poincare Reduction starts with a right (or left) G— invariant La- 
grangian L : TG — > R on the tangent bundle of a Lie group G. Right invariance 
of the Lagrangian may be written as 

L(g(t),g(t)) = L(g(t)h,g(t)h) , for all h G G 

A G— invariant Lagrangian defined on TG possesses a symmetry-reduced Hamil- 
ton's principle defined on the Lie algebra TG/G ~ q. Stationarity of the 
symmetry-reduced Hamilton's principle yields the Euler-Poincare equations 
on the dual Lie algebra q*. For G = Diff , this equation is EPDiff (|2.1|l . 



3 Outline matching & momentum measures 
Problem statement for outline matching: 

Given two collections of curves ci, . . . , cjy and C±, . . . , Cn in fi, find a time- 
dependent diffeomorphic process (t i— > ft) of minimal action (or cost) such that 
ipo = id and y?i(cj) = Ci for i = 1, . . . ,N. The matching problem for the image 
outlines seeks singular momentum solutions which naturally emerge in the 
computation of geodesies. 

3.1 Image outlines as Singular Momentum Solutions of EPDiff 

For example, in the 2D plane, EPDiff has weak singular momentum solu- 
tions that are expressed as |25l23l26j 



where s is a Lagrangian coordinate defined along a set of N curves in the 
plane moving with the flow by the equations x = Q a (i, s) and supported on the 
delta functions in the EPDiff solution (|3.1fl . Thus, the singular momentum so- 
lutions of EPDiff represent evolving "wavefronts" supported on delta functions 
defined along curves Q a (i, s) with arclcngth coordinate s and carrying momen- 
tum P a (t, s) at each point along the curve as specified by l|3.1|> . These solutions 
exist in any dimension and they provide a means of performing CA matching 
for points (landmarks), curves and surfaces, in any combination. 

3.2 Here is the Geometry Leading to the Numerics 

The basic observation that ties everything together in n— dimensions is the fol- 
lowing: 

Theorem (Holm and Marsden, L 23l): EPDiff singular momentum so- 
lutions T*Emb(5, R n ) — ■> g* : (P, Q) — > m define a momentum map. 

It is beyond our scope here to explain either the proof of this theorem or 
the mathematics underlying momentum maps for diffcomorphisms. However, we 
summarize the main results for template matching, as follows: 

— The embedded manifold S is the support set of the P's and Q's. 

— The momentum map is for left action of the diffeomorphisms on S. 

— The whole system is right invariant. 

— Consequently, its momentum map for right action is conserved. 

— These constructions persist for a certain class of numerical schemes. 

— They apply in template matching for every choice of norm. 

3.3 A familiar example of a momentum map 

A momentum map J : T*Q i— > g* is a Hamiltonian for the canonical action 
of a Lie group G on phase space T*Q. It is expressed in terms of the pairing 




(3.1) 



(•,•>: 0* x g R as 



{3,£) = (p, £tq)=:{q*P,Z), 

where (q,p) £ T*Q and the Lie derivative £^q is the infinitesimal generator of 
the action of the Lie algebra element £ £ g on q in the manifold Q. 

The standard example is £%q = £ x g for K 3 xl 3 i— > R 3 , with pairing ( • , • } 
given by scalar product of vectors. The momentum map is then 

J-£=p-(;Xq = qxp-£^J = qxp 

This is angular momentum, the Hamiltonian for phase-space rotations. The 
outlines of images may be parameterized as curves whose dynamics must be 
invariant under rcparametcrizing the arclengths that label those curves. This 
symmetry leads to a conserved momentum map, called the circulation along 
the curves. The analog of this conservation law for fluids is the classical Kelvin 
circulation theorem. 

3.4 EPDifF dynamics informs optimal control for CA 

CA must compare two geometric objects, and thus it is concerned with an op- 
timal control problem. However, the initial value problem for EPDiff also 
has important consequences for CA applications. 

— When matching two geometric structures, the momentum at time t=0 
contains all required information for reconstructing the target 
from the template. This is done via Hamiltonian geodesic flow. 

— Being canonically conjugate, the momentum has exactly the same dimension 
as the matched structures, so there is no redundancy. 

— Right invariance mods out the relabeling motions from the optimal solution. 
This symmetry also yields a conserved momentum map. 

— Besides being one-to-one, the momentum representation is defined on a lin- 
ear space, being dual to the velocity vectors. 

This means one may, for example,: 

• study linear instability of CA processes, 

• take averages and 

• apply statistics to the space of image contours. 

The advantage is the ease of building, sampling and estimating statistical 
models on a linear space. 

3.5 Summary 

We have identified momentum as a key concept in the representation of im- 
age data for CA and discussed analogies with fluid dynamics. The fundamental 
idea transferring from fluid dynamics to CA is the idea of momentum maps 
corresponding to group actions. 



4 Numerical examples of outline matching 



In this section we describe our new technique applying particle-mesh methods to 
the problem of matching outlines. First we describe the approach to calculating 
geodesies in the space of outlines. 

Let Q Q and Q 1 be two embeddings of S 1 in R 2 which represent two shapes, 
each a closed planar curve. We seek a 1-parameter family of embeddings Q(t) : 
S 1 x [0, 1] — > M 2 so that Q(0) = Q and Q(l) matches Q 1 (up to relabeling). 
Q(t) is found by minimizing the constrained norm of its velocity. To find the 
equation for Q we require extremal values of the action 

A = f )-L{u)dt+( I P(s,t)-(Q(s,t)-u(Q(s,t)))dt, L=\\u(t)\\l 



o 2 



i.e. we seek time-series of vector fields u(t) which are minimized in some norm 
subject to the constraint that Q is advected by the flow using the Lagrange 
multipliers P (which we call momentum). The minimizing solutions are 

^ = J P( s ,t)S(x-Q(s,t))ds, (4.1) 

P(s, t) = - P(s, t) ■ Vu(Q(s, t),t), (4.2) 
Q(s,t)=u(Q(s,t),t), (4.3) 

subject to Q(s,0) = Q {s). 

We note that equation (|4.1|) is the momentum map corresponding to the 
cotangent-lift of the action of vector fields u on embedded curves given by 

Q i > u(Q). 

For a suitable test function w, we obtain 

d 8L 
— (w, m) — (Vio, um) + (w, (Vm) t • m) =0, rn — , 
at 5u 

which is the weak form of the EPDiff equation. 

Now one must seek initial momentum P(s,0) which takes shape Qo(s) to 
shape Q 1 (s). To do this, we choose some functional J of the advected shape 
Q(l,s) which is minimized when Q(l, s) matches Q 1 (s). Following [2Hj, we de- 
scribe the curves by singular densities: 

/!=/ fl(s)6(x-Q(l,s)dsdV(x), (4.4) 

77= f r 1 (s)S(x^Q 1 (s)dsdV(x), (4.5) 

and write J = — where || • \\q is a norm for a densities in a reproducing 
kernel Hilbert space with kernel G. This approach means that we do not need 
to force particular points to be matched to each other on the shapes. This last 
problem can be solved by using a gradient algorithm, where the gradient of the 
residual error with respect to P(s,0) is calculated using the adjoint equation 



4.1 Numerical discretization 



We use the Variational Particle-Mesh (VPM) method |30I31| to discretize the 
equations H4.ll4.3f) . as follows: discretize the velocity on an Eulerian grid with 
n g points and approximate \\u\\ there; replace S 1 by representing the shape by 
a finite set of n p Lagrangian particles {Qp}p P = i, and interpolate from the grid 
to the particles using basis functions 

n g n g 

u (Qp) = ^2 u k^k(Qp) > with y^^ fc (cc) = 1, V X. 

k=l k=l 

The action for the continuous time motion on the grid then becomes 




and one can obtain a fully discrete method by discrctizing the action in time. 
For example, we can obtain a first-order method by extrcmizing 

n=l \ \ k ) ) 

The resulting time-stepping method is the (first-order) symplectic Euler-A 
method for the time-continuous Hamiltonian system for the Lagrangian par- 
ticles. In general, the method will always be symplectic since it arises from a 
discrete variational principle (see |32| for a broad introduction to symplectic nu- 
merical methods and their conservation properties). The conservation properties 
of VPM are discussed in (5J . 

We approximate the densities /i and r\ on the grid using the standard particle- 
mesh approach (see [55] 1: 

P 

where Q 1 ^ are the positions of particles on the target shape. This amounts to 
"pixcllating" the singular densities (|4. 414.5(1 on the grid. For a given kernel G, 
we approximate J with 

J = y^,G(x k - xi){fi k - r) k )(m -tn). 

kl 

The discrete adjoint is then applied in computing the inversion for the initial 
conditions for Pp which generate the flow. A numerical example calculated using 
this method is given in figure ^ 
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Fig. 1. Results from a VPM calculation to calculate the minimal path between 
a two simple shapes. On the left, the initial and final shapes are shown, and on 
the right, the deformation of the initial shape into the final shape is depicted 
together with a grid which shows how the flow map deforms the space around 
the shape. We used the H 1 norm for velocity on a 2n x 2-7T periodic domain on 
a 128 x 128 grid, discretized using FFT, and the corresponding kernel was used 
to calculate J. Cubic B-splines were used as basis functions. 
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